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A longstanding conjecture by Belinskii, Khalatnikov, and Lifshitz that the singularity in generic 
gravitational collapse is spacelike, local, and oscillatory is explored analytically and numerically in 
spatially inhomogeneous cosmological spacetimes. With a convenient choice of variables, it can be 
seen analytically how nonlinear terms in Einstein's equations control the approach to the singularity 
and cause oscillatory behavior. The analytic picture requires the drastic assumption that each spatial 
point evolves toward the singularity as an independent spatially homogeneous universe. In every 
case, detailed numerical simulations of the full Einstein evolution equations support this assumption. 
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Introduction. Even though much of the theoretical effort in gravitational physics is focused on either quantum gravity 
or astrophysical sources of gravitational waves, the most fundamental properties of general solutions to Einstein's 
" equations remain unknown. Foremost among these is the nature of the singularity that forms in generic gravitational 



'""^ . collapse. Powerful theorems due to Penrose and Hawking and others |10 1^ predict that some kind of singular 
^ ■ behavior will result. However, the proofs of these theorems require very little of the details of the field equations and 
^ . are thus unable to provide details of the nature of the singularity. Many different types of singular behavior arise in 
' a variety of special cases, but which, if any, of these characterize generic collapse is not known. In fact, it is hard to 
. see that one could hope to answer the question. Einstein's equations are second order nonlinear partial differential 
■ equations (PDE's) in 3 + 1 spacetime dimensions for the metric variables restricted by four constraints and with the 
1^ freedom to impose four coordinate conditions. It is by no means clear that numerical relativity can help to answer 
this question either since a code will "blow up" if any variables become singular. 

The surprising message of this essay then is that progress in understanding singularities in generic gravitational 
collapse can indeed be made. It requires a judicious choice of variables and a synergy between numerical simulations 
which suggest the patterns and analytic predictions which guide the simulations. Our conclusions are based on the 



I , properties of solutions with one or two directions of spatial symmetry. While these are not the most general solutions 
of Einstein's equations, we believe they are sufficiently generic and sufficiently consistent with each other to indicate 
how completely general solutions are likely to behave. 
^ • Approximately thirty years ago, Belinskii, Khalatnikov, and Lifshitz (BKL) proposed, based on a perturbative 
' analysis, that the singularity in generic collapse was spacelike, local, and oscillatory. In their picture, each spatial 
. point evolves toward the singularity as a distinct Mixmaster universe. Their claim has been criticized but never 
' unambiguously proven or refuted. In addition substantial effort is required to follow their arguments. Our recent 
work, however, substantially supports their predictions and suggests a simple scenario: It appears that, despite their 
complexity, Einstein's equations in the approach to the singularity have an attractor which is dominated by nonlinear 
terms acting as effective potentials for the dynamics. These potential terms act primarily on one of the gravitational 
field degrees of freedom and lead to a generic asymptotic solution that supports the BKL description of the generic 
singularity as spacelike, local, and oscillatory. Other information, such as spatial dependence and the behavior of the 
other degrees of freedom, is carried along by the dominant behavior. 

Our method of studying a given family of solutions as they approach the singualrity is essentially the following : 
We first determine the velocity-term dominated (VTD) solutions corresponding to the family. These are solutions 
of the ode's which result from the truncation of Einstein's equations to eliminate all terms containing spatial 
derivatives (plus certain others; see [Q). Note that the VTD solutions are parametrized by spatially-dependent 
constants of integration (the "VTD parameters"). We then substitute these solutions into the Einstein equations or 
into the corresponding Hamiltonian. If, for generic VTD parameters, the terms in the equations or the Hamiltonian 
containing spatial derivatives become exponentially small (hence, dynamically irrelevant) as t — > oo (t = oo is the 
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time of the spacelike singularity), then we predict that the family is asymptotically velocity-term dominated (AVTD) 
in the sense that a generic solution of the Einstein equations asymptotically approaches a VTD solution. If not, 
we monitor which terms might become important (as a result of exponential growth), and we try to predict the 
subsequent behavior. 

We have applied this "Method of Consistent Potentials" (MCP) to various families of cosmological solutions such 
as the Gowdy spacetimes the Gowdy spacetimes with magnetic fields |^^, and U{1) symmetric spacetimes 

[p7|,|6|. In all cases, using a choice of time based on geometric invariants connected with the symmetry, we get a 
clear prediction: either that the spacetime solutions in the family are AVTD, or that they have oscillatory behavior 
characterized by a sequence of AVTD-like epochs punctuated by "bounces" which involve the short-term awakening 
of certain of the "potential" terms. Moreover, numerical simulations of solutions of the full Einstein equations (for 
each family) have in all cases agreed (at generic spatial points) with the MCP predictions. While these numerical 
simulations cannot be followed all the way to the singularity, they can be taken sufficiently far, with sufficient accuracy, 
to verify the simple BKL scenario. 

We shall now discuss some of the details of these studies, starting with a look at spatially homogeneous cosmologies. 

Homogeneous cosmologies. For a spatially homogeneous cosmology with flat spacelike hypersurfaces (Bianchi I), 
the Hamiltonian for Einstein's theory takes the form 

H = -pI+pI+p^_ = 0, (1) 

where pQ, p± are respectively canonically conjugate to Q, the logarithm of the isotropic scale factor, and P±, the 
(logarithmic) anisotropic shears, and we have the Kasner solutions |p^ . These are parametrized by the constants 
v± — p±/\pn\- Asymptotically {Q ~oo), the shear satisfies 

/3± = -v±n (2) 

with v± satisfying as a consequence of the constraint (Q) 

vl + vl = 1. (3) 

For more general homogeneous cosmologies, the Hamiltonian H adds a term V{^l,f3±) which in certain cases jl^] 
takes the form 

T/(r!,/3±)=Cie'^i("^2/3+)_^^^ga.(n+;3++V3/3_)^^^ga3(n+/3+-V3/3„)^ ^ (4) 

where c; and are positive constants. The solutions corresponding to H with this form of the potential are the 
well-studied Mixmaster solutions . To see the oscillatory behavior of Mixmaster solutions from the MCP 

point of view, we note that the Bianchi I Hamiltonian (|^) is exactly the VTD version of the Mixmaster Hamiltonian, 
and the VTD solutions are the Kasner spacetimes. If we now substitute these Kasner solutions into (^, as ^ — oo 
we have 

V -> cie"il"l("^"^"'+) + c2e''^l"l(~^+"'++^''-' -I- c3e"^l"l(-i+''+-^''-' =Vi + V2 + V3. (5) 

A few moments of reflection will show that the restriction (^) implies that any generic values of v± will cause one of 
the Vj to grow and the others to decay. The growing term provides an exponential potential for the dynamics. A 
bounce off the potential changes the sign of the appropriate combination of momenta. What happens then is that a 
pattern of growing and decaying terms continues indefinitely. It has long been known that this picture is borne out 
by the full solutions to Einstein's equations An example showing the evolution of the V's and /3+ vs is given 

in Fig. 1. 

Inhomogeneous cosmologies. We first consider the Gowdy family of spacetimes on x R. These allow spatial 
dependence in one direction with periodic boundary conditions and are described in terms of |]7|, p^ , p0[ P{9, r), Q{9, r), 
essentially the -I- and x polarizations of gravitational waves, and a background A(0, r). This model has the nice feature 
that the dynamical equations for the wave amplitudes P, Q decouple from the constraints which yield A as a function 
of the wave amplitudes. 

Einstein's evolution equations are obtained from the Hamiltonian density ||l6| , ^0|] 

At = ^vTp + e ttq + e F + e ^ 'Q j (b) 
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where the tt's are the conjugate momenta and ' = 9/(90.0 The Hamiltonian constraint is the equation for A (where 
the overdot means d/dr) while the momentum constraint is one for A'. Neglecting spatial derivatives in the evolution 
equations yields, in the limit t ^ oo, the VTD solution 

P = v{0)t, Q^Qo{0), X = -v\0)t, (7) 

with all the momenta dependent only on 9 and v = np / 2tt\ . Substitution in (|^) shows that the only terms in Ti. 
which could fail to be exponentially small are 



v.- ('4) ^ (%) e-, V, . g m e«-.'. (8) 



Vi will grow exponentially if f < 0. V2 will grow if w > 1. Bounces off Vi and V2 eventually drive v into the range 
< w < 1, yielding a consistent AVTD solution at almost every spatial point [|[D. Fig. 2 shows how this process 
works for Vi and V2 and P vs r at a typical spatial point. 

The "magnetic Gowdy" model is obtained by a change of topology and the addition of a magnetic field The 
VTD solution (0) is unchanged. The Hamiltonian (||) picks up an additional term 

V3^A7:xB^e^^+-^/' (9) 

where the constant B measures the strength of the magnetic field and A is now dynamical. Substitution of the VTD 
solution into yields as r ^ 00 

Vs^inlB^e^^-^"^^/^ (10) 

The range < < 1 which allows Vi and V2 to become exponentially small causes V3 to grow. Thus, in the MCP 
picture, no value of v leads to AVTD behavior. Further MCP qualitative study ||2^ indicates that at each spatial 
point, Vi, V2, and V3 alternately awaken, causing bounces characteristic of local Mixmaster behavior. Numerical 
simulations support this prediction with typical behavior shown in Fig. 3. 

It is possible to generalize the plane symmetric Gowdy model on x R by removing one of the Killing vectors. 
This yields vacuum U{1) symmetric models p7|,p|,|6[ . Here, the dynamical degrees of freedom, now called (p and w, no 
longer decouple from the background described by A, z, and x. Einstein evolution equations may be obtained from 
the variation of (,a — d/dx°' with spatial coordinates it, v, and time t) 

H = J J dudv (^-pI + \e''pl + '-pI + \e'-pl - \p\ + 2p^ 

+e-2- f [ dudv{ {e^e'^") ,ab - {e^e^') ,a A,b [{e-^') ,„ - (e'^^) x,u ' 



= J J dudv [n + 2pa) (11) 

where p,^ is canonically conjugate to ^, eab{z, x) is the conformal metric in the u-v plane, and 7i = is the Hamiltonian 
constraint. The choice pA > defines the collapsing direction |6). The VTD solution as t — > 00 is 

z = -WzT, x = xo, A = (2 - p\)t, ip = -v^T, u = Wo, (12) 

where v^, Sq, p\^ v^, lo^ and the remaining momenta are functions of u and v but constants in r. Upon substitution 
of the VTD solution (pj) into (|ll|), it is seen immediately that the Gowdy-like potential terms 



^In the usual treatment of this model P,H], the non-dynamical tta is set equal to i. 
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become negligible if Vz > v^p > 0, as in Gowdy. The very complicated terms in (|11|) with spatial derivatives have a 
simple exponential behavior. Given > 0, the inverse conformal metric e""^ is dominated by e~^^. Examination of 
the terms in (O) shows that all but one of them have the exponential factor 



The remaining term is 



V2 = \e-^^+^e''''uj,a uJb ~ F{x, Vc^)e(-^'A+2-^+4-.)-. (I5) 

Now as T — > 00, the Hamiltonian constraint, with the VTD solution substituted in, becomes 

H^-UPAf + 2vl + 2vl^0. (16) 

This restriction means that > 2vz so that p^, pz > implies the exponential decay of terms with the behavior of 
(p^).pl On the other hand, V2 will decay exponentially only if (p^)^ > (2wz + 4w,^)^ which is inconsistent with jl^). 
Thus V2 must increase exponentially if Vz and are positive. The Hamiltonian constraint restriction means that 
V2 will grow approximately as e'^'"^'^ so that V2 will act as a potential for the ip degree of freedom. A bounce will 
yield < which will cause Vi from (|l|) to grow exponentially. A bounce off Vi will cause the cycle to repeat. All 
other variables will remain consistent with the VTD solution. Surprisingly, this extremely simple picture describes 
the results of numerical simulations of the full U{1) symmetric model Einstein equations (21]. The behavior of Vi, 
V2, and vs T at a typical spatial point is shown in Fig. 4. 

Conclusions. The asymptotic approach to the singularity of a sequence of increasingly complicated cosmological 
spacetimes has been examined. The VTD solution is substituted into the relevant Hamiltonian to obtain a prediction 
as to whether terms not present in the VTD solution become exponentially small. In the homogeneous Mixmaster 
case, we find that no matter what the value of the VTD constants, one of three dominant exponential potential terms 
will always grow. In the vacuum Gowdy case, there are two such potentials at every spatial point and a window, 
< v{9) < 1, consistent with AVTD behavior. Addition of a magnetic field (with a suitable change of topology) 
removes this window and causes the system to oscillate indefinitely as the three potentials grow and decay. Detailed 
study shows this to be truly local Mixmaster behavior |^^. On the other hand, generic U{1) symmetric models have 
oscillations between two rather than three potentials at every spatial point although no value of Vip{u, v) is consistent 
with the AVTD solution. 

We emphasize that, in the inhomogeneous case, this simplified picture follows from the MCP analysis if the time 
dependences of a variety of spatially dependent functions are assumed to be negligible compared to the exponentials 
we have singled out. In every case, this drastic assumption holds true in numerical simulations of the full Einstein 
evolution equations. This means that we could evolve toward the singularity using ODE's at every spatial point rather 
than PDE's. Thus the BKL conjecture that the generic approach to the singularity is local is very strongly supported. 
However, it is not yet clear whether, generically, the further conjecture that our observed oscillations are Mixmaster 
is also correct. 

For the inhomogeneous models, the local nature of the dynamics is plausible. It is easy to show that the horizon 
size goes to zero as r — > 00 in both the Gowdy [|j and U{1) cases. Thus, in the chosen time coordinate, the spatial 
points causally decouple. Information from other sites, carried through the spatial derivatives, ceases to propagate. 

There are several implications of these results. First, the local nature of the approach to the singularity means 
that our choice of cosmological boundary conditions is merely a convenience. The indicated behavior should arise 
asymptotically in any generic collapse leading to a spacelike singularity. Second, as bounces occur at different spatial 
points at different times, the amount of spatial structure in the gravitational waves will increase p^,p|. Third, these 
results are certainly consistent with a generic spacelike singularity (at r = cw) although this was not known before 
these studies to be a property of U{1) symmetric models. 

The simplicity of our picture despite the generality of the one Killing field spacetimes suggests that further gen- 
eralization is possible. Within the models discussed here, one may for example introduce a scalar field of the type 
known to suppress Mixmaster oscillations in the homogeneous case Implementation of the MCP analysis leads 
in all classes of inhomogeneous models studied here to the prediction that the oscillations will be suppressed. These 
predictions have not yet been checked by numerical simulation. 

Finally, the generalization to the zero Killing field case may be within reach. While precise, quantitative numerical 
studies in 2 + 1 and 3+1 dimensions must still come to terms with numerical difficulties caused by the spiky features 



^Numerical simulations show that polarized (7(1) models with V2 = have the predicted AVTD behavior 
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seen in the Gowdy models ||5|,|3|] , the local nature of the collapse means that the qualitative features can be seen even 
at coarse spatial resolutions. 



ACKNOWLEDGEMENTS 



This work was supported in part by NSF Grants PHY9507313 and PHY9722039 to Oakland U., PHY9308177 to 
the U. of Oregon, and PHY9503133 to Yale U. BKB thanks the Institute for Geophysics and Planetary Physics of 
Lawrence Livermore National Laboratory for hospitality. Some numerical simulations were performed at NCSA (U. of 
Illinois). 



Belinskii, V. A., Lifshitz, E. M., and Khalatnikov, I. M., Sov. Phys. C/sp.,13, 745-765, (1971). 

Berger, B. K., Ann. Phys. (N.Y.), 83, 458, (1974). 

Berger, B. K., and Garfinkle, D., Phys. Rev. D, 57, 4767 (1998). 

Berger, B. K., Garfinkle, D., and Strasser, E., Class. Quantum Grav., 14, L29-L36, (1997). 
Berger, B. K., and Moncrief, V., Phys. Rev. D, 48, 4676, (1993). 
Berger, B. K. and Moncrief, V., Phys. Rev. D, bf 57, 15 June 1998. 
Gowdy, R. H., Phys. Rev. Lett, 27, 826, (1971). 

Grubisic, B., and Moncrief, V., Phys. Rev. D, 47, 2371-2382, (1993). 
Halpern, P., J. Gen. Rel. Grav., 19, 73-94, (1987). 

Hawking, S. W., and Penrose, R., Proc. Roy. Soc. Lond. A, 314, 529-548, (1970). 

Hobill, D., Burd, A., and Coley, A., eds.. Deterministic Chaos in General Relativity, (Plenum, New York, 1994). 

Isenberg, J. A., and Moncrief, V., Ann. Phys. (N.Y.), 199, 84, (1990). 

Kasner, E., Trans. Am. Math. Soc., 27, 155-162, (1925). 

Kirillov, A. A., and Kochnev, A. A., JETP Lett, 46, 435-438, (1987). 

Misner, C. W., Phys. Rev. Lett, 22, 1071-1074, (1969). 

Moncrief, V., Ann. Phys. (N.Y.), 132, 87-107, (1981). 

Moncrief, V., Ann. Phys. (N.Y.), 167, 118, (1986). 

Ryan, Jr., M. P., and Shepley, L. C, Homogeneous Relativistic Cosmologies, (Princeton University, Princeton, 1975). 
Wald, R. M., General Relativity, (University of Chicago Press, Chicago, 1984). 
Weaver, M., Isenberg, J., and Berger, B. K., Phys. Rev. Lett 80, 2984 (19 98). 



Berger, B. K. and Moncrief, V., submitted to Phys. Rev. D, gr-qc/9804085 



5 



-40 




100 200 300 400 500 



- 



•12 



l°9ioV 



-24 



600 



FIG. 1. A typical Mixmaster trajectory (vacuum Bianchi IX) from a numerical simulation of Einstein's equations. The 
potentials Vj are defined in (0). 
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FIG. 2. Evolution at a representative spatial point in a numerical simulation of Einstein's equations for a typical vacuum 
Gowdy model. The potentials Vi and V2 are defined in dsl). 
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FIG. 3. Evolution at a representative spatial point in a numerical simulation of Einstein's equations for a typical magnetic 
Gowdy model. The potentials Vi and V2 are defined in (0) while V3 is defined in (S). 




FIG. 4. Evolution at a representative spatial point in a numerical simulation of Einstein's equations for a typical generic 
U{1) symmetric model. The potentials are defined in ( |l^ ) and ([l^). 
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